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We investigate the low-energy properties of the Holstein polaron through calculation of the g-dependent 
phonon spectral function using an improved exact-diagonalization technique, defined over a variational Hilbert 
space. We perform a comprehensive study of the low-energy excitations of the polaron. Beside the energy range, 
where the additional phonon excitation is unbound, we observe separate coherent peaks which correspond to 
bound and antibound states of a polaron and additional phonon quanta. These novel states can be observed 
for intermediate and strong electron-phonon coupling strengths, as well as below and above the unbound one- 
phonon excitation spectrum. A detailed investigation of their properties is presented. We find good agreement 
between the phonon spectral function obtained from the first-order strong-coupling perturbation theory and 
numerical results. 
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I. INTRODUCTION 

Electron-phonon coupling represents one of the fundamen- 
tal mechanisms that determines thermodynamic as well as 
transport properties of solids. Polaron formation is a phe- 
nomenon where a single charge carrier changes its state by 
absorbing or emitting phonon quanta. Its importance has been 
identified in many novel materials 1 , such as high-temperature 
superconductors 2-4 , colossal magnetoresistance materials 5 , 
molecular crystals and fullerenes^. In the case where the 
electron couples to an optical branch of lattice vibrations, a 
widely accepted approximation is used where the electron- 
phonon (e-ph) coupling strength and the phonon energy are 
considered momentum independent. This leads to a Holstein 
molecular crystal model, which is one of the most fundamen- 
tal models in solid state physics. Many analytical and numer- 
ical methods have provided a fair understanding of its ground 
state properties, i.e., the formation of a polaroid—. Its behav- 
ior, particularly in the intermediate-coupling (IC) and strong- 
coupling (SC) regime, strongly depends on the strength of 
the phonon energy ojq. In the antiadiabatic regime where the 
Lang-Firsov transformation 13 in the SC limit provides a com- 
prehensive picture of the polaron, most approximations usu- 
ally manage to calculate the system properties qualitatively 
correctly, while the treatment of adiabatic quantum phonons 
poses a more challenging task and is a matter of the recent 
scientific interest^—. 

Above the polaron ground state energy E°, several features 
of the excited states can be observed. In particular, we are in- 
terested in the low-energy excited states which do not belong 
to the one-phonon continuum of an unbound phonon starting 
at E° + ojq. A while ago, the emergence of a state below 
the one-phonon continuum has been reporte d 17 ' 18 . In addi- 
tion, a coherent state above the continuum has been found us- 
ing dynamical coherent potential approximation 19 . In the last 
decades, the investigation of these states spread to all param- 
eter regimes of the Holstein model. In the SC adiabatic ap- 
proximation 2 ^ 2 ^, several states were observed below the con- 



tinuum. When moving to non-zero phonon energies, one such 
state was found in the SC approximation 22 and three states in 
the IC regime using exact-diagonalization technique 23 . Re- 
cently, an investigation of properties of the coherent states 
below the one-phonon continuum was performed by vari- 
ational exact diagonalization 24 ' 25 and momentum-averaging 
(MA) approach 2 ^. 

Regarding dynamical properties of the Holstein polaron, 
the analysis rarely focused on the states outside the contin- 
uum. Spectral weight was assigned to these states in the 
calculation of one-electron spectral function using dynami- 
cal mean-field theory (DMFT) 2 ^ and MA approach 28 , optical 
conductivity using DMFT 30 , and in the quantum Monte Carlo 
study of the one-electron spectral function in the Rashba- 
Pekar model 3 -^ Less emphasis has been devoted to the phonon 
spectral function. The latter was studied in Ref. 32 , where au- 
thors numerically and analytically calculated the gr-dependent 
phonon spectral function for various phonon energies and for 
all relevant e-ph coupling regimes. In the low-lying excita- 
tion spectrum they find a mirror peak, which was supposed to 
represent a state lying an energy ujq above the polaron peak. 
They did not focus on the description of the excited states be- 
low the one-phonon continuum. On the other hand, the author 
of Ref. 33 studied lattice correlation functions of few excited 
states below the continuum in the IC regime, representing the 
spectral weight of the phonon spectral function. Calculating 
these correlations, the author showed that these states contain 
a non-vanishing spectral weight and should contribute to a for- 
mation of coherent excited bands. Nevertheless, the features 
of the spectrum at and above ojq were not studied. Taking this 
into account, we argue that a comprehensive study of the low- 
lying features of the Holstein polaron phonon spectral func- 
tion is still missing. In particular, the emergence of the states 
below and above the one-phonon continuum, denoted here as 
the novel states, needs to be characterized. 

In this paper we study the low-energy spectrum of the Hol- 
stein polaron model, reflected through calculation of the q- 
dependent phonon spectral function. Using an efficient exact- 
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diagonalization technique defined over a variational Hilbert 
space 24 , we are able to investigate a polaron band as well as 
a non-dispersive continuum band starting at ljq above the po- 
laron ground state. This band represents the states where the 
additional phonon excitation is not bound to the polaron. In 
addition, the emergence of novel states with the energy below 
and above the one-phonon continuum is studied. These states 
are denoted in the text as bound and antibound states, respec- 
tively. Calculation of the static correlation function in these 
states shows that the weight of the extra phonon excitations 
decreases exponentially with the distance from the polaron. 
For different values of the model parameters, a particular em- 
phasis is given to the onset of the antibound state, not studied 
before in the literature. The emergence of one bound and an- 
tibound peak in the phonon spectral function is studied within 
the framework of the first-order strong-coupling perturbation 
theory. For a certain range of the model parameters, this ana- 
lytical calculation provides results that are in good agreement 
with the numerical solution. 

The paper is organized as follows. In Sec. II. we introduce 
the model and the numerical method. In Sec. III. we show 
the numerical results, while in Sec. IV. we compare the re- 
sults with the first-order strong-coupling perturbation theory. 
A summary is given in Sec. V. 

II. MODEL AND NUMERICAL METHOD 

We start by writing the one-dimensional spinless Holstein 
model as 

H = -t 22(c\cj + H.c.) +g m (at + q f ) + u }^ qto,-, 

(1) 

where c\ and a\ are electron and phonon creation operators 
at site i, respectively, and rt,; = cjcj electron density at site i 
(with rii = 0, 1). The total number of electrons in the system 
is equal to 1. ojq denotes a dispersionless optical phonon en- 
ergy and t nearest-neighbor hopping amplitude (we set t = 1 
in Sec. I. -III.). The e-ph coupling strength is denoted as g and 
will be in the following replaced by a dimensionless coupling 
9 = 9/uo- 

The phonon Green function, that corresponds to the 
displacement-displacement correlation function, is defined as 

D{q,tt-t2) = -i(i> \f[x q {h)xl{t2)]\Tp°), (2) 

where x q — a q + a ! _ q and T the time ordering operator. Ap- 
plying the Fourier transform and calculating the spectral rep- 
resentation B(q,uj) = —^lmD R (q,ui), where w > and 
D R (q,Lu) the retarded Green function, one gets the corre- 
sponding phonon spectral function 

B q (u) = ]T \{^\x^°)\ 2 5(u - [E% - (3) 

n 

with E and tp° denoting the polaron ground state energy and 
wave function, respectively, while the sum runs over all ex- 
cited states. The ground state wave function if) has momen- 
tum k = 0, and q may be nonzero. 



We used an improved numerical method, originally intro- 
duced in Ref. 24 , which constructs the variational Hilbert space 
(VHS) starting from the single-electron Bloch state 0^10) on 
infinite lattice. The VHS is then generated by applying the 
off-diagonal terms of Hamiltonian ([U to the starting state 

{\^ N r M) )} = (H kin +H^) Nh cim, (4) 

where iJkin and H g corresponds to the first and the second 
term of the Hamiltonian in Eq. \T\ respectively. The set of 
basis states is determined by parameters Nh and M, where 
Nh — 1 is the maximum distance between the electron and 
the phonon quanta and Nh * M is the maximum number of 
phonon quanta contained in the Hilbert space. The param- 
eter M > 1 first introduced in Ref. 34 , is chosen to ensure 
good convergence for the strong e-ph coupling in the adiabatic 
regime where the low-energy states contain multiple phonon 
excitations. On the other hand, large values of Nh provide 
a fair description of the polaron states also in the weak cou- 
pling regime, where the spacial extent of a lattice deforma- 
tion around the electron is large. We have used Nh > 8 and 
M > 4 that lead to converged results in all parameter regimes. 
Note that for any Hilbert space in our calculations the number 
of sites is infinite. The convergence towards thermodynamic 
limit is then achieved when the number of phonon states for 
a given electron location is sufficiently increased. Once the 
VHS is generated, the Holstein Hamiltonian is diagonalized 
using standard Lanczos procedure, where translational sym- 
metry is explicitly taken into account. 

III. RESULTS 

To start the calculation of the phonon spectral function, de- 
fined in Eq. [3] we first focus on the lowest-energy peak that 
represents the polaron ground state. Its spectral weight is 
proportional to the electron density 1 /N and increases when 
approaching the SC limit as g 2 . In the adiabatic regime 
cjo *C 1 the cross-over from large to small polaron occurs 
for A = e p /2t > 1, where e p = g 2 cuo is the polaron energy at 
t = and A a dimensionless e-ph coupling. In Fig. 02a) the 
g-dependent phonon spectral function is plotted for luq = 0.2 
and A = 1.05. At this e-ph coupling strength the bandwidth 
of the polaron W is reduced to W = (Wo = Wo/ 10, where 
Wo is the polaron bandwidth in the weak-coupling (WC) limit 
(Wo — ojo if ojo < 4 and Wo = 4 otherwise) and ( a di- 
mensionless parameter. When ujo is increased away from the 
adiabatic regime, the cross-over to SC regime is less abrupt. 
In Fig. [l~fb) it can be seen for luq — 1.0 and A = 1.5 that the 
polaron dispersion is still substantial and ( = 0.31. 

We now turn to the investigation of the low-lying excited 
states, representing the main focus of the paper. Naively, the 
low-lying excited states should consist of a polaron at a par- 
ticular fc-point and an additional unbound phonon excitation 
with momentum Q, leading to an arbitrary total momentum 
q = k + Q of the composite excited state. This implies that in 
the thermodynamic limit of the phonon spectral function, the 
unbound excited states would have a non-vanishing spectral 
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Figure 1: (Color online) (a) g-dependent phonon spectral function 
B q {co) for q = 0, ...tt and cj = 0.2, A = g 2 /2tuj = 1.05. Only 
the polaron peak and one-phonon excitation spectrum is shown. We 
have used an artificial damping of e = 0.003. The Roman numerals 
I., II., III. denote the three bound states, (b) B q (uf) foruo = 1-0 and 
A = 1.5. Here only the lowest bound state (Roman numeral I.) can 
be observed. An artificial damping of e = 0.01 was used. The Ro- 
man numeral IV. denotes the antibound peak, where the one-phonon 
excitation has an energy higher than ujq + W. The non-dispersive 
peak (denoted by R) corresponding to one-phonon excitation repelled 
from the polaron, is located slightly above ljo + W. This peak is de- 
scribed in more detail in Fig. [2] and Fig. [4] and further in the text. 
The grey areas correspond in (a) and (b) to a continuum of states in 
the interval [ujq , cjo + W] ■ We used the Hilbert space obtained for 
N h = 8, M = 11 that lead to N st = 5.9 x 10 5 . 



weight in the interval [ujq, uiq + W]. In Fig.Q] this continuum 
of states starting at uj is shaded and its width is given by the 
bandwidth W of the polaron peak. Due to the numerical re- 
strictions of the variational Hilbert space, the calculations do 
not show a considerable spectral weight throughout the whole 
continuum band, even though its approximate width is well 
reproduced. When increasing the Hilbert space, the density 
of states inside the continuum band grows. The effects of the 
finite-size Hilbert space will be discussed later. 

Beside the unbound one-phonon excitations, it is well 
known that in the IC regime of the Holstein model a bound 
state of a polaron and additional phonon excitations appears, 
where the additional phonon excitations are bound to the po- 
laron and the energy of such composite state is less than loq 
above the ground state. This phenomenon is sometimes de- 
noted as softening. It has been shown 24,25 that the average 
value of phonon quanta in the bound state N^ h = a\ai) 
can not be related to the corresponding value of the ground 
state as N^ h = N® h + 1. The latter holds true only for the 
unbound one-phonon excitations. In general, there are two 
bound states with different symmetry of the wave functions. 
In Fig. [Ha) forwo = 0.2 and A = 1.05, the lowest bound state 
denoted by the Roman numeral I. has a symmetric wave func- 
tion at q = and an antisymmetric wave function at q = ir, 
while the second bound state (the Roman numeral II.) is sym- 
metric at q — tt. Their spectral weight in B q (ui) is approx- 
imately proportional to electron density 1/N and vanishes at 
q = state 33 . In the case of A = 1.05, an additional ex- 



cited state (the Roman numeral III.) appears below the phonon 
threshold u . This state, however, has a non- vanishing spec- 
tral weight also at q = 0. The emergence of the bound peaks 
in B q (u>) is best resolved in the IC regime as shown in Fig.Q] 
since with the increasing of e-ph coupling towards the SC 
regime, their positions approach the phonon threshold ujq. On 
the other hand, when the phonon energy is increased (the case 
for ujq = 1.0 is shown in Fig.[TJb)) the positions of the bound 
peaks shift as well closer to the bare phonon energy. In this 
case only the lowest bound state (the Roman numeral I.) can 
be resolved from the spectrum. 

While the non-zero spectral weight in B q (uj) for the 
peaks below the phonon threshold was suggested in previous 
works 33 , we find a novel coherent peak located above the con- 
tinuum band of unbound one-phonon excitations. This peak 
is denoted by the Roman numeral IV. in Fig. |TJb). In analogy 
to the bound states, this peak represents a state where the po- 
laron and additional phonon excitations are antibound, i. e., its 
energy is above the upper edge of the continuum band of the 
unbound one-phonon excitations, E q > ojq+W. Note that for 
ujq = 0.2, as seen in Fig. [TJa) no antibound peak can be ob- 
served as its spectral weight is vanishingly small. One of the 
first indications that the energy spectrum of the Holstein po- 
laron can contain a state above the one-phonon continuum for 
the non-adiabatic values of the phonon energy, was reported 
many years ago in the calculation based on dynamical coher- 
ent potential approximation 19 . Recently, coherent states above 
the one-phonon continuum were observed in the study of one- 
electron spectral function using DMFT 29 and MA approach 28 . 

The antibound peak is more pronounced for the higher val- 
ues of luq. In Fig. |2{a), (b) and Fig. |2jc) we plot the q- 
dependent phonon spectral function for ojq = 2.0 and 6.0, 
respectively. As the strength of the dimensionless e-ph cou- 
pling g 2 approaches the SC limit, the spectral weight of the 
antibound peak increases. The energy of its state relative to 
the ground state E q is the highest at q = tt, however its dis- 
persion is not proportional to the polaron dispersion. When 
the momentum is decreased, the energy decreases, but always 
remains above the continuum band. Its spectral weight van- 
ishes at the q = point, where the symmetry of the antibound 
wave function is antisymmetric. On the other hand, the wave 
function is symmetric at q = it. Note that in all cases under 
investigation when ljq > 1, we use the dimensionless e-ph 
coupling g = g/cuQ. In the antiadiabatic limit ljq ^> 1, the 
condition g ^ 1 determines the cross-over from large to small 
polaron. 

The phonon spectral function calculated in Ref. 32 using the 
kernel polynomial method can be compared to our results. We 
find a similarity between the antibound peak as denoted in our 
paper and the mirror peak as defined in the latter reference. 
Even though the authors of Ref. 32 did not examine the struc- 
ture of this peak, we may speculate that these two peaks rep- 
resent the same state. In the second part of this section, we 
shall further investigate the properties of the antibound state 
through calculation of the e-ph correlation function (Fig. [U 
and its evolution when the e-ph coupling increases (Fig. [5] and 

In addition to the bound and antibound peaks, described 
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Figure 2: (Color online) (a), (b) g-dependent phonon spectral func- 
tion B q (Lu] I for g 2 = 1.0 (where ( = 0.62) and g 2 = 2.0 (C = 0.29), 
respectively, and u)o ~ 2.0. (c) g-dependent phonon spectral func- 
tion B q (uj) for g 2 — 1.5 (£ = 0.22) and ujo = 6.0. The grey areas 
correspond to a continuum of states in the interval [ujq , + W] ■ An 
artificial damping of e = 0.02 was used in all cases. The Hilbert 
space used in the calculation was obtained for Nh = 11, M = 8. 
The inset of (c) shows B q (uj) around wo > 6.0 at g = 7r for dif- 
ferent system sizes. The lowest peak (red curve) corresponds to 
a Hilbert space (N h , M) = (8,11), N st = 5.9 x 10 5 . As the 
weight of the peak increases and its position shifts to lower ener- 
gies, Hilbert space is changed to (10, 7), (12, 7) and (14, 4), which 
implies N st = 1.8 x 10 6 , 2.5 x 10 7 and 1.8 x 10 7 , respectively. 



above, another feature of the phonon spectral function can be 
observed in Fig. Q] and Fig. [2ja) and (b), denoted by the letter 
R. This is a non-dispersive peak located slightly above ojq + W 
with the spectral weight comparable to the peaks due to bound 
and antibound phonon excitations. However, due to the non- 
dispersive nature of the peak, its origin should be different 
from the latter peaks. Even though it seems from Fig. Q] and 
Fig.|2]that its spectral weight is not sensitive to the value of u>q, 
it suddenly vanishes as uo > 4.0, as demonstrated in Fig.[2J c ) 
for wo = 6.0. We will show later in the calculation of the e-ph 
correlation function (Eq. |5j that this state is formed by a po- 
laron and an additional phonon excitation, where the latter is 
repelled from the polaron in real space. We shall argue as well 
that the position of this peak approaches ujq + W from above 
as the size of the Hilbert space increases and in the thermody- 
namic limit merges with the continuum. In contrast, spectral 
weight and positions of peaks I. to IV. are well converged in 
the thermodynamic limit. 

The energy of the polaron ground state calculated by our 
numerical method is variational in the thermodynamic limit as 
Nh —> oo and can be calculated to a great accuracy 35 . This is 
also reflected in the calculation of the polaron energy in differ- 
ent q-points, where its peak position is not sensitive to differ- 



ent sizes of Hilbert space (not shown in the paper). For the un- 
bound one -phonon excitations, where the additional phonon 
excitation is not correlated to the polaron position and can 
be located anywhere on the lattice, the energy convergence 
is slower. The dominant contribution of the unbound phonon 
excitations in the phonon spectral function comes from the 
peak at luq above the polaron ground state. In our calculations 
the position of this peak is located at ujq + e for a finite Nh, 
where e decreases when Nh is increased. For a fixed system, 
e is the largest in the WC limit. In the current work, our main 
focus is the behavior of the phonon spectral function in the 
IC regime. Generation of a sufficient number of phonon ex- 
citations in the vicinity of the charge carrier in this regime, 
starts to play a more crucial role on the system properties than 
the maximal distance between polaron and the extra phonon 
excitation. While the first condition is well controlled by the 
value of the parameter M in the generation of VHS, the lat- 
ter is determined by the parameter Nh- We use this subtle 
interplay between the values of these two parameters generat- 
ing the Hilbert space to achieve the convergence of both un- 
bound, bound and antibound low-energy states in the phonon 
spectral function. The inset of Fig.[2|d) show the position of 
the wo-peak for q — ir, g 2 = 1.5 and cjo = 6.0. As Nh in- 
creases from Nh = 8 to 14, the peak indeed approaches the 
bare phonon energy and gains the spectral weight. We have 
thus shown that for various sizes of the Hilbert space included 
in our calculation, the wo-peak is well reproduced. 

It is crucial to show that not only the wo-peak, but the bound 
and antibound peaks as well as the width of the continuum 
band converge well within the framework of our numerical 
method. This can be seen from Fig.[3j where the low-energy 
excitation spectrum of B q {uj) is plotted for ojq = 2.0 and 
g 2 = 3.0. Results are obtained for different system sizes 
Nh = 10, 12 and 14 while the parameter M is fixed to 4. 
There are two important features in this figure. The first one 
concerns the continuum of the unbound phonon excitations in 
the interval [ujq,ujo + W]. The density of the states in this 
interval depends on the maximal number of allowed phonon 
quanta in our calculation as well as on the maximal allowed 
distance between the phonon quantum and the polaron. As the 
parameter Nh is increased, the unbound states become denser, 
leading in the thermodynamic limit to the continuum of states 
with the non-zero spectral weight. On the other hand, the 
lowest bound peak and the antibound peak are located out- 
side the continuum. Since in the case of the bound and an- 
tibound states the additional phonon quanta are attached to 
the polaron, their position and weight in the phonon spectral 
function should be less sensitive to the finite Hilbert space. In- 
deed, the three nearly overlapping curves for different Hilbert 
spaces in Fig. [3] indicate that these peaks are well converged 
states, clearly separated from the continuum of states. This is 
not surprising since in this case, as we shall see in Fig. [4] the 
amplitude of the extra phonon excitations decrease exponen- 
tially with the distance from the polaron and can be described 
efficiently within the VHS. 

In order to investigate the structure of the observed peaks 
of the B q (uj) in more detail, we calculate the static electron- 
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Figure 3: (Color online) B q (uj) for g 2 = 3.0 and ojq — 2.0 
(£ = 0.12). The polaron peak is not shown and the repulsive peak 
is removed from the figure. The grey area corresponds to the contin- 
uum band of unbound states in the interval [cjo , ojo+W^] . An artificial 
damping of e = 0.007 was used. The three curves were obtained us- 
ing the following parameters (Nh, M): (10, 4) - red, (12, 4) - yellow 
and (14, 4) - black, that lead to the following number of basis states: 
7V st = 2.2 x 10 5 , 2.0 x 10 B and 1.8 x 10 7 . The peaks below and 
above the grey area, where the curves are nearly overlapping, cor- 
respond to the lowest bound state (green) and to the antibound state 
(blue), respectively. 
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Figure 4: (Color online) Electron-phonon correlation number 7, (i — 
j) for the ground state and various excited states at ojq = 2.0. (a) 
Polaron ground state at q — and q = tt for g =1.0. (b) An 
excited state with momentum q — it and g 2 — 1.0, which lies an 
energy ljo above the state in the inset of (a), (c) The lowest bound 
state at q — and q — tt/2 for g 2 — 3.0. (d) Antibound and 
repulsive state at q — it and g 2 = 1.0. Note that the integral for 
the repulsive state over the two humps marked by arrows, is equal 
to 1.0. The insets show the correlation functions on a logarithmic 
scale. The Hilbert space used in these calculations was obtained for 
(N h , M) = (12, 7), N st = 2.5 x 10 7 . 



phonon correlation function 

^(i-j) = {^\4c i a]a j \^), (5) 

which represents the distribution of the number of excited 
phonon quanta in the vicinity of the electron. The strength 
of the numerical method is reflected in the fact that Jg(i — j) 
can be calculated for (i — j) > 10 in any state q as well as for 
any low-lying excited state n obtained by the Lanczos diago- 
nalization. 

Using 7™ (i—j) of Eq. (O, different properties of the excited 
states can be analyzed. For example, the unbound state can be 
distinguished from the bound and antibound state. In Fig. [4] 
we plot 7™(i — j) at uq = 2.0. We start with the polaron state 
at g 2 = 1.0, Fig. SJa). The ground state at the q = point 
has the minimum number of phonon quanta, located in the 
vicinity of the electron, while the number of phonon quanta 
increases when q approaches the Brillouin-zone boundary. To 
show the spatial correlation of the electron and phonon posi- 
tion, the exponential decay of 7°(« — j)\ q=7 v is shown in the 
inset of Fig. Ufa). The one-phonon excitation of the latter state 
is shown in Fig.HJb). As seen from the inset of the figure, this 
state does not exhibit an exponential decay, i.e., the additional 
phonon excitation is unbound. 

It is interesting to compare Fig.Hfb) with (c) and (d), where 
we plot 7™(i — j) for the first bound (at g 2 — 3.0) and the 



antibound state (at g 2 = 1.0), respectively. Note that in Fig. [3] 
the finite size investigation was performed for u>q — 2.0 at 
g 2 = 3.0 due to the poor resolution of the bound peak for the 
lower values of the e-ph coupling (compare with Fig. [2 a) and 
(b)). Nevertheless, the antibound peak detaches at the q = it 
point from the continuum already at g 2 ~ 0.77. In the inset 
°f Fig-Elc) and (d), 7™(i — j) is shown on a logarithmic scale 
for q = it/ 2 (bound state) and q = tt (antibound state), re- 
spectively. For both cases the exponential decay confirms that 
in such composite excited states, the extra phonon excitations 
are bound to the polaron. 

Other important information can be obtained from 7™(i — 
j). As already mentioned in the discussion of Fig. Q] and [2] 
there is a non-dispersive peak in B q {ui) denoted by R with 
a considerable spectral weight located above the continuum 
of unbound phonon excitations. When increasing the size of 
the Hilbert space (not shown in the paper), its spectral weight 
does not scale to zero and remains a well defined quantity. 
In Fig. |Sd), the structure of this state can be resolved for 
(Nh, M) = (12, 7). Its striking features are two humps, well 
separated from the polaron peak. The integral over each of the 
humps gives 0.5 for any momentum, which implies that an ad- 
ditional phonon excitation is repelled from the polaron. Con- 
sequently we call the non-dispersive peak obtained in Fig. Q] 
and |2j the repulsive peak. Due to the repulsion between the 
polaron and an extra phonon excitation, its energy will always 
remain above luq + W in the finite Hilbert space calculations. 
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Thus the energy gap between this peak and the continuum 
band is an artifact of the method and scales to zero in the 
thermodynamic limit. Consequently, the repulsive peak would 
merge with the upper edge of the continuum band of the un- 
bound states. 



B (co) 

q 



cg =6.0 




Figure 5: (Color online) B q (uS) for different values of e-ph coupling 
strength, (a) q = and (b) q — n, both at ujq = 2.0. (c) q = and 
(d) q = 7r, both at uio = 6.0. For the latter phonon energy, the po- 
laron peak is not shown. An artificial damping of e = 0.02 was used. 
The Hilbert space was obtained for parameters (Nh, M) = (8, 11). 
The blue peak, denoted by the Roman numeral IV., represents the 
antibound peak, while the red peak dashed by horizontal lines and 
denoted by R represents the repulsive peak. 

By now, our aim has been to identify and introduce all the 
distinct features of the phonon spectral function and to prove 
their robustness to the finite-size Hilbert space. We would 
now like to get some general insight into the emergence of 
these novel states for different values of the phonon energy 
and the e-ph coupling strength. 

In Fig. [5] we plot the phonon spectral function for a fixed 
q and different values of g 2 . We focus here mainly on the 
antibound and repulsive peak. In Fig. |5J a ) an d (b), B q (ui) is 
shown for q = and q = ir, respectively, and ojo = 2.0. At 
these values of parameters, the particular g-points do not show 
any spectral weight of the bound states. On the other hand, 
the repulsive peak denoted by R located slightly above luq + 
W, is clearly visible and gains its spectral weight when the 
e-ph coupling is increased. The peaks in the continuum band 
of the unbound one-phonon excitations display a maximum 
at ojo, while their distribution and spectral weights decrease 
towards the upper edge luq + W. Nevertheless, due to the 
emergence of the repulsive state in our calculations, we would 
anticipate an increase of the spectral weight at the upper edge 
of the continuum band. When ujq is increased to 6.0, as seen 
in Fig. |3c) and (d), the repulsive state and correspondingly 



the nonmonotonous decay of the unbound excitation spectrum 
disappears. 

The antibound peak observed in the same spectra for q = n 
denoted by the Roman numeral IV. is detached from the con- 
tinuum of states and gains as well its spectral weight when the 
e-ph coupling is increased. In addition, it is more pronounced 
for higher values of ujq when calculated for the same values 
of the dimensionless e-ph coupling g 2 . This should be com- 
pared to Fig. [Tfb) where B q (uj) is plotted at ujq = 1.0 and 
A = 1.5 (i.e., g 2 — 3.0). In this figure the antibound peak 
can be hardly observed. Note that there is the repulsive peak 
in Fig. Q2b) located above the well converged antibound peak. 
As the system size Nh is further increased, the repulsive peak 
moves to lower energies, as indicated before. 
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Figure 6: (Color online) Energy positions of different states in the 
phonon spectral function above the ujq threshold vs. the strength 
of the e-ph coupling. Dashed line marks the state at uio above the 
polaron ground state, while the dotted line marks the state at 2u)q. 
Bound states are not shown in the figure. The black curve denotes the 
wo + W state and the blue curve the antibound state at q = n with 
the energy E^. Vertical dashed-dotted line marks the lowest value of 
g 2 , at which the antibound peak can be observed at a particular ujq. 
The Hilbert space was obtained for parameters (Nh, M) = (8, 8) 
that lead to N st = 2.0 x 10 5 . 

In general, there is no theoretical argument against emer- 
gence of additional coherent states above the ljq + W thresh- 
old in the energy spectrum of the Holstein polaron. In Sec. 
IV. we calculate the one-phonon excitation spectrum of the 
polaron in the SC antiadiabatic limit, as described within the 
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Lang-Firsov picture. The first order perturbation theory pre- 
dicts only one state above the uj + W continuum in this limit, 
consistent with our numerical results. As seen from Fig.|5]and 
|6l this state, denoted throughout the paper as the antibound 
state, persist up to IC regime of e-ph coupling and moderate 
values of u>q. However, in the adiabatic limit ujq -C 1 where 
multiple excited coherent states exist below the ojq continuum, 
it is likely that several coherent states emerge above the con- 
tinuum as well. In terms of the phonon spectral function as 
shown in Fig. QI a X the latter states have vanishingly small 
spectral weight and are not discussed in the paper. 

When addressing the emergence of the antibound state in 
B q {uj), a question arises at which values of e-ph coupling the 
corresponding peak can be resolved from the spectrum. Since 
in the last part of this section (Fig.|7| we show that the onset of 
this peak exhibit a non-monotonous behavior when the model 
parameters are changed, we need to get some more insight 
into the properties of the antibound state. 

In Fig.[6]the energy positions relative to the polaron ground 
state energy E° for different states above the luq threshold 
are plotted vs. the strength of the e-ph coupling. The gray 
area corresponds to the continuum band with the upper edge 
ujq + W and the blue curve, i.e., the full line emerging out of 
the continuum, represents the antibound state with the energy 
E% at q = 7r. A particular emphasis is given to the onset of 
the antibound state. A vertical dashed-dotted line marks the 
lowest value of g 2 , for which this state can be observed. It 
emerges out of the one-phonon continuum if luq > 2.0 (see 
Fig. Hfa)-(c)), or out of higher-energy continuum if ujq < 2.0 
(see Fig. Hd)). 

The above analysis shows that the antibound peak is not 
part of the continuum band and its energy is always higher 
than ujq above the top of the polaron band. Thus we propose a 
slightly different interpretation of this peak than in Refi^. In 
our picture the peak corresponds to the state where the addi- 
tional phonon excitations are attached to the polaron. This is 
qualitatively different from the states in the continuum band, 
where the extra phonon excitation is not bound to the polaron. 
In Sec. IV. we will derive similar results from the calculation 
of B q (u>) within the first-order strong-coupling perturbation 
theory. 

To conclude the discussion we present a diagram in Fig. [7] 
For a particular value of usq, it shows the minimal strength of 
the e-ph coupling that the bound or antibound state would de- 
tach from the band of the unbound one-phonon excitations. 
For the antibound state, the minimum e-ph coupling is in the 
state q = ir, while for the lowest bound state the minimum 
appears at q = 0. Note that for the bound peak, its spectral 
weight is zero at q = 0. Consequently, the corresponding 
curve in Fig. |7]represents just its lower boundary. It is in good 
agreement with Re f. 24-26 , where both analytical approxima- 
tions as well as numerical calculations are compared. As in 
Fig- different regimes of wo and g are included in the dia- 
gram, not all of the parameters are of physical interest. For a 
fixed ljq, a state at ( = W/Wq = 10 _1 is represented by a 
dotted line and a state at£ = 10~ 2 bya solid line. Above the 
latter value, the polaron is self-trapped and can be in the non- 
adiabatic regime well described by the Holstein Hamiltonian 



in the atomic limit t = 0. 




Figure 7: (Color online) The minimum strength of the e-ph coupling 
for a fixed value of ujq, at which the bound or antibound state de- 
tach from the band of the unbound one-phonon excitations. The an- 
tibound state is plotted for q = -n and the bound state for q = 0. 
Note that the e-ph coupling energy g was introduced in Eq. [TJ The 
dashed line represents the case when g — ujq (i.e., the result from the 
first-order sept for both bound and antibound state), while the dotted 
and solid line correspond to the cases when the polaron bandwidth is 
reduced to W = (Wo- The inset shows the same results in terms of 
the dimensionless adiabatic parameter A. 

As already noticed in Fig. [6] where the vertical dashed- 
dotted lines mark the onset of the antibound peak in B q (uj), 
the corresponding curve in Fig. [7] shows a non-monotonous 
behavior. It can be divided in two parts connected to each 
other at ujq 2.0. While in the adiabatic regime the spectral 
weight of the antibound state is vanishingly small (see Fig. [T} 
and its emergence begins at the e-ph coupling values compa- 
rable to that of the bound states, the most interesting behavior 
starts at wo > 2.0. In particular, for 2.2 < wo < 3.0, the onset 
of the antibound peak occurs at g 2 < 0.5, i.e., for the values 
of e-ph coupling in the IC regime. The dashed curve in Fig. [7] 
represents the case when g = 1.0, which is as well the onset 
of the bound and antibound state derived from the first-order 
strong-coupling perturbation theory (sept). Comparing this to 
the numerical results in the non-adiabatic regime, we see that 
for a fixed ujq, the emergence of the bound state is shifted to 
higher values of the e-ph coupling while the emergence of the 
antibound state occurs always at g < 1.0. 

IV. STRONG-COUPLING PERTURBATION THEORY 

The phonon spectral function can be calculated analytically 
in some limiting cases. Starting from the SC limit t = where 
the Hamiltonian of Eq. Q] describes a displaced harmonic os- 
cillator on the site of the electron, we derive the expression 
for B q (uj) in the first order perturbation of hopping t. We fol- 
low the procedure initiated in Ref. 24,36 . Our aim is to show 
that both bound and antibound state can be described already 
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by the 1st order sept and their spectral weight in B q (ui) com 
pared to the numerical result. 

Applying the canonical transformation H = < 

S = -gT, j ( a j- a 'j) n j 

V, where the unperturbed part of H is 



s He~ s with 
one gets the Hamiltonian H = Hq + 



H 



w 

3 



E'< 



(6) 



and V is considered a perturbation 

V = -ter 92 ^(c]c j+1 e- 9{a h a hi) e g ^- a ^ + H.c). 

3 

(7) 

We have introduced the new basis where a harmonic oscillator 
has a shifted origin on the electron site j. We denote such state 
as \gj). The ground state wave function at arbitrary g-point, 
which corresponds to the energy E q — — e p — 2te 9 cosq, 
can be written in translationally invariant form as \ijj q ) = 
l/V^V^Zj e%q ^\9j)i where N is the number of lattice sites. 
The bosonic operators dj and at in Eq. [6] and [7] act relative to 
Ifli) if j — ii an d relative to an unshifted oscillator if j ^ i. 

Once the ground state wave function is known, the polaron 
spectral weight in B q (uj) can be obtained from the matrix el- 
ement 



«I4 



N 



which leads to the polaron peak 

B° q {uj) = ^-6 (u - 2te- g \\ - cos q) 



(8) 



(9) 



The weight of the peak is proportional to the e-ph coupling 
parameter g 2 . In Fig.[8j B q (uj) is shown at wo/t = 2.0 and 
g 2 = 2.5, calculated (a) from the first-order sept and (b) nu- 
merically. The dispersion and the bandwidth of the polaron 
peak show good agreement since ( yields 0.16 in (a) and 0.18 
in(b). 

The main focus of the calculation is the low-energy exci- 
tation spectrum. There are N degenerate wave functions that 
consist of the translationally invariant electron state and an ad- 
ditional phonon quantum at a fixed distance d away from the 
electron 



i 



E< 



\9j 



i j+d\y3i- 



(10) 



The definition of the bosonic operators is the same as in Eq.|6] 
and [7] To obtain the energies of the excited states one needs 
to calculate the matrix elements V^ d , = {i'WV]^,) . For 
g > 1, there are N — 2 solutions lying in the energy interval 
[E^,EY\, where E± = —e p +ujQ±2te~ 9 , and two additional 
states lying below and above the interval. These are the bound 
and antibound state, respectively. The emergence of these two 
states at g = 1 is plotted by the dashed curve in Fig. [7] For 
g > 1, the bound state first detaches from the continuum at 



q = and the antibound state at q — n. When g is further 
increased as shown in Fig.[8ja) for g 2 = 2.5, these two states 
detach from the continuum for all values of q. In the g — > 
oo limit, both bound and antibound state approach the energy 
-e p +uj . 
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Figure 8: (Color online) B q (co) for uio/t — 2.0 and g 2 = 2.5. (a) 
B q (ui) — Bq(ui) + Bq(uj), calculated from the first-order sept for 
iV = 51 sites (Eq.|9]and Eq. \\2\. (b) Numerical results, obtained 
using the parameters (Nh, M) = (8, 11). The bound peak (the Ro- 
man numeral I.) and the antibound peak (the Roman numeral IV.) 
can be observed in both cases. The repulsive peak (denoted by R) is 
obtained only in (b). An artificial damping of e = 0.01 was used. 



The spectral weight of any one extra phonon excitation in 
,(lu) can be obtained from the matrix element 

(tfMW°) = -7fiT, b *b> i9d > (H) 



where \t/j q ) — J2dbd(q)\?Pd) is the corresponding eigenstate 
of the matrix V 1 . Note that (^|a_ g |^°) = for any state. 
The total spectrum of the one extra phonon excitations is than 



E & ™ 



5(u-(El-E )). (12) 



The comparison of the phonon spectral function between 
the first-order sept and numerical results is shown in Fig. [8] 
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for g 2 — 2.5 and uio/t = 2.0. While the repulsive peak does 
not appear in the first-order sept approximation, the bound and 
antibound peak can be resolved from both spectra at very sim- 
ilar energy positions and comparable spectral weights. The 
wave function of the bound state is symmetric at q = and 
antisymmetric at q = it, while the opposite is true for the 
antibound state. This indicates that the essential physics of 
these two novel states is well captured by the first-order sept, 
where the calculation of the one extra phonon excitations, as 
determined by the matrix elements V} d , is analogous to the 
ID tight-binding problem, modified around the position of 
the electron 24 . In fact, if one assumes a regular tight-binding 
model with a non-zero onsite energy eo at only one site, its 
analytical solution yields an additional localized state below 
(for eo < 0) or above (for eo > 0) the continuous band of 
itinerant solutions. 



V. SUMMARY 

In conclusion, the main focus of the paper has been to cal- 
culate the low-energy features of the phonon spectral func- 
tion for the phonon energies ujq ranging from adiabatic to an- 
tiadiabatic regime. Using an improved exact diagonalization 
method defined over a variational Hilbert space, we were able 
to get converged results in a wide range of the e-ph coupling 
strengths with special emphasis on the properties in the IC 
regime. In particular, we were interested in the low-energy 
excitation spectrum of the polaron, where the most spectral 
weight originates from the states where the polaron and the 
additional phonon excitation are unbound. These states are 
limited to the energy interval with the width equal to the po- 
laron bandwidth W. Even though in our calculations we use 
a finite Hilbert space, the results are consistent with the ex- 
istence of a continuum band of unbound excited states in the 
thermodynamic limit. 

In addition, we have identified the emergence of novel 
states in the spectrum of the low-energy excitations for dif- 
ferent regimes of model parameters. Below the continuum 
band we have found three bound states of the polaron and 
additional phonon excitations, with the energy less than uj 



above the polaron ground state. On the other hand, we have 
found one state above the one-phonon continuum denoted as 
the antibound state which consists of the polaron and addi- 
tional phonon excitations with the energy higher than ujo + W 
above the polaron ground state. In both cases we showed that 
the spatial correlation between the polaron and extra phonon 
excitations decrease exponentially with the distance from the 
polaron. Consequently, the numerical method can calculate 
their properties to a sufficient accuracy. Beside the bound and 
antibound states, we identified the state with a considerable 
spectral weight where the additional phonon excitation is re- 
pelled from the polaron. Due to the finite Hilbert space in 
our calculations, the energy of such repulsive state is always 
higher than that of the unbound states. However, we antic- 
ipate in the thermodynamic limit that this peak merges with 
the upper edge of the continuum band. 

The bound and antibound states develop in the IC regime 
of the e-ph coupling and persist up to the SC limit. We have 
performed a detailed investigation of the properties of the an- 
tibound state, especially the onset of the corresponding peak 
in the phonon spectral function for a wide range of the model 
parameters, which lead to an alternative interpretation of this 
peak regarding its former explanation. The antibound peak 
can be observed for the non-adiabatic values of loq, whereas 
the adiabatic limit of the states above the uq + W threshold is 
not analyzed in the paper. Our results can be compared to the 
analytical solution of the first-order sept, valid for the antiadi- 
abatic values of the phonon energies. In this approximation, 
the bound and antibound state can be observed for g > luq. 
Nevertheless, the accurate numerical investigation yields even 
a lower value of the e-ph coupling for the onset of the anti- 
bound state (Fig. [7]), indicating that the antibound peak can 
be already observed before the cross-over to the small polaron 
regime. 
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